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: Abstract 

Pythia version 6 represents a merger of the Pythia 5, Jetset 7 and SPythia programs, 
with many improvements. It can be used to generate high-energy-physics 'events', i.e. sets 
of outgoing particles produced in the interactions between two incoming particles. The 
objective is to provide as accurate as possible a representation of event properties in a 
wide range of reactions. The underlying physics is not understood well enough to give 
an exact description; the programs therefore contain a combination of analytical results 
and various models. The emphasis in this article is on new aspects, but a few words of 
general introduction are included. Further documentation is available on the web. 



CPC subject index: 11.2 

PACS codes: 13.60.-r, 13.65. +i, 13.85.-t, 12.15.-y, 12.38.-t, 12.60.-i 
Keywords: event generators, multiparticle production 



New Version Summary 



Title of program: Pythia 
Version number: 6.154 

Catalogue identifier: — 

Distribution format: uuencoded compressed tar file 

References to most recent previous versions: Computer Physics Communications 82 
(1994) 74 and Computer Physics Communications 101 (1997) 232, respectively 

Catalogue identifiers of most recent previous versions: ACTU 

Authors of most recent previous versions: T. Sjostrand, and S. Mrenna, respectively 
Does the new version supersede the previous versions?: yes 

Computers for which the new program is designed and others on which it has been tested: 
Computer: DELL Precision 210 and any other machine with a Fortran 77 compiler 

Installation: Lund University 

Operating system under which the new program has been tested: Red Hat Linux 6.2 

Programming language used: Fortran 77; is also fully compatible with Fortran 90, i.e. 
does not make use of any obsolescent features of the Fortran 90 standard 

Memory required to execute with typical data: about 800 kwords 

No. of bits in word: 32 (double precision real uses two words) 

No. of processors used: 1 

Has the code been vectorized or parallelized?: no 
No. of bytes in distributed program: about 1.8 Mb 

Keywords: QCD, standard model, beyond standard model, hard scattering, e''"e~ annihi- 
lation, leptoproduction, photoproduction, hadronic processes, high-p^ scattering, prompt 
photons, gauge bosons, Higgs physics, parton distribution functions, jet production, par- 
ton showers, fragmentation, hadronization, beam remnants, multiple interactions, particle 
decays, event measures 

Nature of physical problem: high-energy collisions between elementary particles normally 
give rise to complex final states, with large multiplicities of hadrons, leptons, neutrinos and 
photons. The relation between these final states and the underlying physics description is 
not a simple one, for two main reasons. Firstly, wc do not even in principle have a complete 
understanding of the physics. Secondly, any analytical approach is made intractable by 
the large multiplicities. 

Method of solution: complete events are generated by Monte Carlo methods. The com- 
plexity is mastered by a subdivision of the full problem into a set of simpler separate tasks. 
All main aspects of the events are simulated, such as hard-process selection, initial- and 
final-state radiation, beam remnants, fragmentation, decays, and so on. Therefore events 
should be directly comparable with experimentally observable ones. The programs can 
be used to extract physics from comparisons with existing data, or to study physics at 
future experiments. 

Restrictions on the complexity of the problem: depends on the problem studied 
Typical running time: 10-1000 events per second, depending on process studied 
Unusual features of the program: none 
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1 Introduction 



The Pythia and Jetset programs [|T| are frequently used for event generation in high- 
energy physics. The emphasis is on multiparticle production in colhsions between elemen- 
tary particles. This in particular means hard interactions in e"'"e~, pp and ep colliders, 
although also other applications are envisaged. The programs can be used both to compare 
with existing data, for physics studies or detector corrections, and to explore possibili- 
ties at present or future experiments. The programs are intended to generate complete 
events, in as much detail as experimentally observable ones, within the bounds of our 
current understanding of the underlying physics. The quantum mechanical variability 
between events in nature is here replaced by Monte Carlo methods, to obtain 'correctly' 
both the average behaviour and the amount of fluctuations. Many of the components of 
the programs represent original research, in the sense that models have been developed 
and implemented for a number of aspects not covered by standard theory. 

Although originally conceived separately, the Pythia |^ and Jetset § programs 
today are so often used together that they have here been joined under the Pythia header. 
To this has been added the code of SPythia 0|, an extension of Pythia that also covers 
the generation of supersymmetric processes. The current article is not intended to give 
a complete survey of all the program elements or all the physics — we refer to the long 
manual and further documentation on the Pythia web page 

http : / /www . thep . lu . se/ ~torb j orn/Pythia . html 
for this. Instead the objective is to give a survey of new features since the latest official 
publications, with some minimal amount of background material to tie the story together. 
Many of the advances are related to physics studies, which are further described in separate 
articles. Others are of a more technical nature, or have been of too limited a scope to 
result in individual publications (so far). 

Some of the main topics are: 

• An improved simulation of supersymmetric physics, with several new processes. 

• Many new processes of beyond-the-standard-model physics, in areas such as techni- 
color and doubly-charged Higgses. 

• An expanded description of QCD processes in virtual-photon interactions, combined 
with a new machinery for the flux of virtual photons from leptons. 

• Initial-state parton showers are matched to the next-to-leading order matrix ele- 
ments for gauge boson production. 

• Final-state parton showers are matched to a number of different first-order matrix 
elements for gluon emission, including full mass dependence. 

• The hadronization description of low-mass strings has been improved, with conse- 
quences especially for heavy-flavour production. 

• An alternative baryon production model has been introduced. 

• Colour rearrangement is included as a new option, and several alternative Bose- 
Einstein descriptions are added. 

Many further examples will be given. In the process, the total size of the program code 
has almost doubled in the six years since the previous main publication. 

The report is subdivided so that the physics news are highlighted in section 2 and the 
programming ones (plus a few more physics ones) in section 3. Section 4 contains some 
concluding remarks and an outlook. 
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2 Physics News 



For the description of a typical high-energy event, an event generator should contain a 
simulation of several physics aspects. If we try to follow the evolution of an event in some 
semblance of a time order, one may arrange these aspects as follows: 

1. Initially two beam particles are coming in towards each other. Normally each par- 
ticle is characterized by a set of parton distributions, which defines the partonic 
substructure in terms of flavour composition and energy sharing. 

2. One shower initiator parton from each beam starts off a sequence of branchings, 
such as q — > qg, which build up an initial-state shower. 

3. One incoming parton from each of the two showers enters the hard process, where 
then a number of outgoing partons are produced, usually two. It is the nature of 
this process that determines the main characteristics of the event. 

4. The hard process may produce a set of short-lived resonances, like the Z°/W^ gauge 
bosons, whose decay to normal partons has to be considered in close association with 
the hard process itself. 

5. The outgoing partons may branch as well, to build up final-state showers. 

6. In addition to the hard process considered above, further semihard interactions may 
occur between the other partons of two incoming hadrons. 

7. When a shower initiator is taken out of a beam particle, a beam remnant is left 
behind. This remnant may have an internal structure, and a net colour charge that 
relates it to the rest of the final state. 

8. The QCD confinement mechanism ensures that the outgoing quarks and gluons are 
not observable, but instead fragment to colour neutral hadrons. 

9. Normally the fragmentation mechanism can be seen as occurring in a set of separate 
colour singlet subsystems, but interconnection effects such as colour rearrangement 
or Bose-Einstein may complicate the picture. 

10. Many of the produced hadrons are unstable and decay further. 
In the following subsections, we will survey updates of the above aspects, not in the 
same order as given here, but rather in the order in which they appear in the program 
execution, i.e. starting with the hard process. 

2.1 Physics Subprocesses 
2.1.1 Process Classifications 

Pythia contains a rich selection of physics scenarios, with well above 200 different sub- 
processes, see Tables H ^ and |^. The process number space has tended to become a bit 
busy, so processes are not always numbered logically. Some processes are closely related 
variants of the same basic process (e.g. the production of a neutralino pair in processes 
216-225), others are alternative formulations (e.g. TPtP h° has to be convoluted with 
the flux of Z°'s around fermions, and thus is an approximation to fjfj — > fjfjh°; when 
the latter was implemented the former was still kept). A process may also have hidden 
further layers of processes (e.g. can be produced in top decays, whichever way top is 
produced). 

One classification of the subprocesses is according to the physics scenario. The follow- 
ing major groups may be distinguished: 

• Hard QCD processes, i.e. leading to jet production. 
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Table 1: Subprocesses, part 1: standard model, according to the subprocess numbering 
of PYTHlA.'f denotes a fermion (quark or lepton), 'Q' a heavy quark and 'F' a heavy 
fermion. 



No. Subprocess 



Hard QCD processes: 

11 

12 fjf j — i> ffcffc 

13 idi ^ gg 
28 fig ^ fig 
53 gg ffeffc 
68 gg — > gg 



Soft QCD processes: 

91 elastic scattering 

92 single diffraction {XB) 

93 single diffraction {AX) 

94 double diffraction 

95 low-px production 



Open heavy flavour: 
(also fourth generation) 

81 fif j QfcQfe 

82 gg^QfcQfc 

83 Qfcfz 

84 g7^QfcQfc 
85 



Closed heavy flavour: 

86 gg^J/V'g 

87 gg xocg 

88 gg xicg 

89 gg X2cg 

104 gg ^ xoc 

105 gg-> X2c 

106 gg J/V^7 

107 g7 J/V'g 

108 77 J/^7 



No. Subprocess 



W/Z 

1 

2 
22 
23 
25 
15 
16 
30 
31 
19 
20 
35 
36 
69 
70 



production: 



fA- 
f.f, - 
fA- 
f.f, - 
fA- 
f f — 

fj, - 

fA- 
fj, - 

77 ^ 
7W± 



7VZ0 

zoz° 

w+w- 

. gW± 
fiZO 

7ZO 

. 7W± 

w+w- 



Prompt photons: 



14 

18 fif 

29 
114 
115 



^ g7 
, 77 

f.7 

77 
g7 



Deep inelastic scatt. 

10 fA^fA 

99 7A^/^ 



Photon-induced: 

33 fi7 fig 

34 fi7^f^7 
54 g7 ffcffc 
58 77 f Jfc 



No. Subprocess 



131 fi7T^fig 

132 fi7£^f.g 

133 fi7T fi7 

134 fi7£^f*7 

135 g7^^fA 

136 g7£^fA 

137 7*7* ^fA 

138 7T7L^f*fi 

139 7£7T^f^fi 

140 7£7£^fA 
80 qi7 qfcvr^ 



Light SM Higgs: 



3 


fA- 




24 


fA- 


-.Z^hO 


26 


fj, - 


W±hO 


102 


gg - 


.hO 


103 


77 - 


^ho 


110 


fA- 


^7hO 


121 


gg - 


^ QfeQ..h° 


122 




- Q;^.Q,h° 


123 


ff - 


f.f,hO 


124 


ff - 


ffcf/h" 



Heavy SM Higgs: 



5 


Z^ZO ^ 


hO 


8 


w+w- 


^ho 


71 


Z^ZO^ 


Z^ZO 


72 


Z^ZO-^ 


W+Wl 


73 


z^w^ 


-z^w^ 


76 


W+Wl 


-z^zo 


77 







Soft QCD processes, such as diffractive and elastic scattering, and minimum-bias 
events. Hidden in this class is also process 96, which is used internally for the 
merging of soft and hard physics, and for the generation of multiple interactions. 
Heavy-flavour production, both open and hidden (i.e. as bound states like the J/ip). 
Hadronization of open heavy flavour will be discussed in section |2.4.1|. Some new 
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Table 2: Subprocesses, part 2: beyond the standard model non-SUSY, with notation as 
above. 



No. Subprocess 



BSM 

151 

152 

153 

171 

172 

173 

174 

181 

182 

156 

157 

158 

176 

177 

178 

179 

186 

187 



Neutral Higgses: 
iiU ^ HO 
gg^HO 
77 ^ H° 

ff 

gg Q/cQ^H" 

fif. - AO 
gg^ AO 
77 AO 
UU ^ ZOAO 



f^f.-HO 



fif,AO 



f f 



4f/A0 



gg ^ QfcO^AO 



Charged Higgs: 
143 Ufj H+ 

161 ug-^m' 



Higgs pairs: 

297 Ufj H±hO 

298 Ufj H^RO 

299 fA ^ AOhO 

300 fA ^ AOro 

301 UU H+H- 



No. Subprocess 



New gauge 

141 fA^ 

142 iifj - 
144 fj, - 



bosons: 
7/Z0/Z'° 
■ W'+ 
• R 



Technicolor 
149 gg ^ 



191 
192 
193 
194 
195 
361 
362 
363 
364 
365 
366 
367 
368 
370 
371 
372 
373 
374 
375 
376 
377 



f f 
f f 

f f 
f f 
ff 

if 
f f 
f f 

f f 
f f 
f f 
f f 
f f 
f f 



Vtc 
Pt°c 

Pte 



TTtcTTtc 

ZO^t 



TT. 



tc 



w 



W^ZO 
W^vrO 

0^± 



ZOtt. 



tc 



W^ttO 



No. Subprocess 



Compositeness: 



146 
147 
148 
167 
168 
169 
165 
166 



e7 — 
dg - 
ug - 

qiqj 
qiqj 
qiqi 
f^f.(- 



e" 
d* 
u* 
d*qk 

e±e*^ 
. 77Z0) 
^ W±) - 



ffcffc 



Doubly-charged Higgs: 
341 



342 
343 
344 
345 
346 
347 
348 
349 
350 
351 
352 



f.f, - 



. H^±e^ 
' H^^e^ 

^ H^P^ 

. H^r^ 
H^HZ- 
H^^H^" 



ffcfzH 
f.f,H 



Leptoquarks: 

145 qiij Lq 

162 qg-^£LQ 

163 gg LoL, 

164 q^q^ 



LqLq 



processes have been added for closed heavy flavour, but we remind that data here 
are yet not fully understood, and have given rise to models extending on the more 
conventional Pythia treatment 0. 

W/Z production. A first-order process such as fjfj gW^ is now quite accurately 
modeled by the initial-state shower acting on fjfj — > W^, see section |2.2.1| , but the 
former can still be useful for a dedicated study of the high-p^ tail. 
Prompt-photon production. 

Photon-induced processes, including Deep Inelastic Scattering. A completely new 
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Table 3: Subprocesses, part 3: SUSY, with notation as above. A trailing + on a final 
state indicates that the charge-conjugated one is included as well. 



No. Subprocess 



SUSY: 
201 
202 
203 
204 
205 
206 
207 
208 
209 
210 
211 
212 
213 
214 
216 
217 
218 
219 
220 
221 
222 
223 
224 
225 
226 
227 
228 
229 



U^ 
Q^ 
U^ 
U^ 
kU 

U^ 
^^^^ 
ff 
f f 
ff 
f f 



f.f. 
fj^ 
f.f^ 

f.f. 

f.f. 
U^ 

U^ 



XiXi 

X2X2 
X3X3 

X4X4 

X1X2 
X1X3 

XiXi 

X2X3 
X2XA 

X3X4 

XiXi 

X2X2 
XIX2 

XiXi 



No. Subprocess 



230 
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
246 
247 
248 
249 
250 
251 
252 
253 
254 
256 
258 
259 
261 
262 



f f 
f f 

f f 
ff 
ff 
f f 

fj, 
f.f^ 

f.f^ 
ff 
f f 
f f 



Ug- 
Ug- 
Ug- 
Ug- 
Ug- 

Ug- 

Ug- 
Ug- 
Ug- 
Ug- 



X2X1 

Xsxt 
XiXi 
Xixt 

X2X2 
X3X2 

XiXt 
gXi 

gX2 
gX3 
gXA 

gXi 

gxt 



(liLXl 
^iRXl 
^iLX2 
^iRX2 
^iLX3 
^iRX3 

qii?X4 

CLjLXi 
CLjLxt 

CLiLg 
CLiRg 

>tit^ 

• t2t2 



No. Subprocess 



263 
264 
265 
271 
272 
273 
274 
275 
276 
277 
278 
279 
280 
281 
282 
283 
284 
285 
286 
287 
288 
289 
290 
291 
292 
293 
294 
295 
296 



f,;f,: 



f f 

f f 

hi-] 

f f 

hi- J 

f f 

hi-] 

f f 

hl-j 

fj, 
f^f. 
f,;f. 



bqi 
hqi 
bqi 
bqi 
bqi 
bq, 

qiq, 
q^q* 



bb 
bb 

bb 

bg 
bg 
bb 



^tit^+ 
till 

' qiLqjL 

^ (liR^ljR 
' ^iLqjR+ 

* qiLq^L 

^ (iiR(l*jR 

^ qiL%*R+ 
' qjiq^L 

' (ijR^ljR 
(liR^iR 

* biqii 
^ bzqii?, 

biQiij, + b2qiL 
' biq*L 
' h2q*R 

^ biqjij + b2q*L 

bib* 
-^bab; 
bib* 

h2hl 
bibi 
b2b2 
bib2 
big 

b2g 
bib;+ 



machinery for 7*p and 7*7* physics has been constructed here, see section 2.1.3 



Standard model Higgs production, where the Higgs is reasonably light and narrow, 
and can therefore still be considered as a resonance. 

Gauge boson scattering processes, such as WlWl WlWl (L = longitudinal), 
when the standard model Higgs is so heavy and broad that resonant and non- 
resonant contributions have to be considered together. 

Non-standard Higgs particle production, within the framework of a two-Higgs- 
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doublet scenario with three neutral (h", H° and A°) and two charged (H^) Higgs 
states. Normally associated with SuSY (see below), but does not have to be. The 
Higgs pair production processes were previously hidden in process 141, but are now 
included explicitly. 

• Production of new gauge bosons, such as a Z', W and R (a horizontal boson, 
coupling between generations). 

• Technicolor production, as an alternative scenario to the standard picture of elec- 
troweak symmetry breaking by a fundamental Higgs. Processes 149, 191, 192 and 
193 may be considered obsolete, since the other processes now include the decays to 
the allowed final states of the pl^/u^^/pf^ bosons, also including interferences with 
7/Z°/W^. The default scenario is based on [^. 

• Compositeness is a possibility not only in the Higgs sector, but may also apply to 
fermions, e.g. giving d* and u* production. At energies below the threshold for new 
particle production, contact interactions may still modify the standard behaviour; 
this is implemented not only for processes 165 and 166, but also for 11, 12 and 20. 

• Left-right symmetric models give rise to doubly charged Higgs states, in fact one 
set belonging to the left and one to the right SU(2) gauge group. Decays involve 
right-handed W's and neutrinos. The existing scenario is based on |^. 

• Leptoquark (Lq) production is encountered in some beyond-the-standard-model sce- 
narios. 

• Supersymmetry (Susy) is probably the favourite scenario for physics beyond the 
standard model. A rich set of processes are allowed, even if one obeys i?-parity 
conservation. The supersymmetric machinery and process selection is inherited 
from SPythia [Q, however with many improvements in the event generation chain. 
Relative to the SPythia process repertoire, the main new additions is sbottom 
production, where a classification by mass eigenstates is necessary and many Feyn- 
man graphs are related to the possibility to have incoming b quarks. Many different 
Susy scenarios have been proposed, and the program is flexible enough to allow 
input from several of these, in addition to the ones provided internally. 

Obviously the list is far from exhaustive; it is a major problem to keep up to date with 
all the new physics scenarios and signals that are proposed and have to be studied. 
One example of another physics area that has attracted much attention recently is the 
possibility of extra dimensions on 'macroscopic' scales. Also, a general-purpose program 
can not be optimized for all kinds of processes. If a generator for some kind of partonic 
configurations is already available, outside of Pythia, there exists the possibility to feed 
this in for subsequent treatment of showers and hadronization. 



2.1.2 Parton Distributions 



For cross section calculations, the hard partonic cross section has to be convoluted with 
the parton distributions of the incoming beam particles. The current default is GRV 94L 
for protons and SaS ID for real and virtual photons 0. Some further parameterizations 
are available in Pythia, such as the recent CTEQ 5 proton ones |JlO|], and a much richer 
repertoire if the Pdflib library |11| is linked. 
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2.1.3 Photon Physics 

Since before, a model for the interactions of real photons is available, i.e. for 7p and 77 



events [|12[- This has now been improved and extended also to include virtual photons, i.e. 



7*p and 7*7* events |T^. It is especially geared towards the transition region of rather 
small photon virtualities ^ 10 GeV^, where the physics picture is rather complex, while 
it may be overkill for large Q^, where the picture again simplifies. 

Photon interactions are complicated since the photon wave function contains so many 
components, each with its own interactions. To first approximation, it may be subdivided 
into a direct and a resolved part. (In higher orders, the two parts can mix, so one has to 
provide sensible physical separations between the two.) In the former the photon acts as a 
pointlike particle, while in the latter it fluctuates into hadronic states. These fluctuations 
are of 0{aem), and so correspond to a small fraction of the photon wave function, but 
this is compensated by the bigger cross sections allowed in strong-interaction processes. 
For real photons therefore the resolved processes dominate the total cross section, while 
the pointlike ones take over for virtual photons. 

The fluctuations 7 — > qq {—>■ 7) can be characterized by the transverse momentum k± 
of the quarks, or alternatively by some mass scale m ~ 2k ±, with a spectrum of fluctu- 
ations oc dk'^/kj_. The low-fcj_ part cannot be calculated perturbatively, but is instead 
parameterized by experimentally determined couplings to the lowest-lying vector mesons, 
V = p°, u;", (j)^ and J/ip, an ansatz called VMD for Vector Meson Dominance. Parton dis- 
tributions are defined with a unit momentum sum rule within a fluctuation , giving rise 
to total hadronic cross sections, jet activity, multiple interactions and beam remnants as 
in hadronic interactions. States at larger k± are called GVMD or Generalized VMD, and 
their contributions to the parton distribution of the photon are called anomalous. Given 
a dividing line ko ~ 0.5 GeV to VMD states, the parton distributions are perturbatively 
calculable. The total cross section of a state is not, however, since this involves aspects 
of soft physics and eikonalization of jet rates. Therefore an ansatz is chosen where the 
total cross section scales like ky/k^, where the adjustable parameter ky ~ rnp/2 for light 
quarks. The spectrum of states is taken to extend over a range ko < k± < ki, where ki is 
identified with the p±mm{s) defined in eq. below. There is some arbitrariness in that 
choice, and for jet rate calculations also contributions to the parton distributions from 
above this region are included. 

If the photon is virtual, it has a reduced probability to fluctuate into a vector meson 
state, and this state has a reduced interaction probability. This can be modeled by a 
traditional dipole factor {my/lmy + Q'^))'^ for a photon of virtuality Q^, where my — ^ 2k± 
for a GVMD state. Putting it all together, the cross section of the GVMD sector then 
scales like 

1^1 dkl kl ( Akl V . , 



hi kl kl \4ki + gv ' 

A real direct photon in a 7p collision can interact with the parton content of the 
proton: 7q qg and 7g — > qq. The p± in this collision is taken to exceed fci, in order to 
avoid double- counting with the interactions of the GVMD states. For a virtual photon the 
DIS (deeply inelastic scattering) process 7*q — >• q is also possible, but by gauge invariance 
its cross section must vanish in the limit — > 0. At large Q^, the direct processes 
can instead be considered as the 0{as) correction to the lowest-order DIS process. The 
DIS 7*p cross section is here proportional to the structure function F2{x,Q'^) with the 
Bjorken x = Q'^/{Q'^ + W^). Since normal parton distribution parameterizations are 
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frozen below some Qq scale and therefore do not obey the gauge invariance condition, an 
ad hoc factor {Q^/iQ"^ + '^p))^ is introduced for the conversion from the parameterized 
F2{x,Q'^) = Seqg(x,(5^) to a ctdis- In order to avoid double-counting between DIS and 
direct events, we decide to introduce a requirement p± > max(A;i, Q) on direct events. In 
the remaining DIS ones, thus p± < Q. The DIS rate should be reduced accordingly, by 
a Sudakov form factor giving the probability not to have an interaction above scale Q, 
which can be approximated by exp(— (T2iI.ect/'^DIs)• 
Note that the dependence of the DIS and direct processes is implemented in the 
matrix element expressions. This is different from VMD/GVMD, where dipole factors are 
used to reduce the assumed flux of partons inside a virtual photon relative to those of a 
real one, but the matrix elements contain no parton virtuality dependence. 

After some further minor corrections for double-counting, we arrive at a picture of 
hadronic 7*p events as being composed of four main components: VMD, GVMD, direct 
and DIS. Most of these in their turn have a complicated internal structure, as we have 
seen. The 7*7* collision between two inequivalent photons contains 13 components: four 
when the VMD and GVMD states interact with each other ('double-resolved'), eight with 
a direct or DIS photon interaction on a VMD or GVMD state on either side ('single- 
resolved', including the traditional DIS), and one where two direct photons interact by 
the process 7*7* qq ('direct', not to be confused with the direct process of 7*p). 

Several further aspects can be added to the above machinery. The impact of resolved 
longitudinal photons is unknown, except that it has to vanish in the limit ^ 0, and 
can be approximated by some Q^-dependent enhancement of the normal transverse one. 
For a complete description of ep events or e"'"e~ two-photon ones, a convolution with the 
X- and Q^-dependent flux of virtual photons inside an electron is also now provided. 



2.1.4 Supersymmetry 

Pythia simulates the Minimal Supersymmetric Standard Model (MSSM), based on an 
effective Lagrangian of softly-broken SuSY with parameters defined at the weak scale, 
which is typically between mz and 1 TeV. The MSSM particle spectrum is minimal in 
the sense that it includes only the partners of all Standard Model particles (presently 
without massive neutrinos), a two-Higgs doublet — one Higgs Hu coupling only to up- 
type fermions and one coupling only to down-type fermions — and partners, and the 
gravitino. Once the parameters of the softly-broken SuSY Lagrangian are specified, the 
interactions are fixed, and the sparticle masses can be calculated |0 . 



The masses of the scalar partners to fermions, sfermions, depend on soft scalar masses, 
trilinear couplings, the Higgsino mass fi, and tan (3, the ratio of Higgs vacuum expectation 
values (Hu)/(Hd). The masses of the fermion partners to the gauge and Higgs bosons, 
the neutralinos and charginos, depend on soft gaugino masses, fi, and tan /5. Finally, the 
properties of the Higgs scalar sector is calculated from the input pseudoscalar Higgs boson 
mass rriA, tan/?, fi, trilinear couplings and the sparticle properties in an effective potential 



approach ||T5|. Of course, these calculations also depend on SM parameters (mt,mz,as, 
etc.). Any modifications to these quantities from virtual MSSM effects are not taken into 
account. In principle, the sparticle masses also acquire loop corrections that depend on 
all MSSM masses. 

i?-parity conservation is assumed (at least on the time and distance scale of a typ- 
ical collider experiment), and only lowest order, sparticle pair production processes are 
included. Only those processes with e~^e~ , fi~^fi~ , or quark and gluon initial states are 



9 



simulated. Likewise, only i?-parity conserving decays are allowed, so that one sparticle 
is stable, either the lightest neutralino, the gravitino, or a sneutrino. SuSY decays of the 
top quark are included, but all other SM particle decays are unaltered. 

Various improvements to the simulation are being implemented in stages. Some of 
these can have a significant impact on the collider phenomenology. Among these are: 
the generalization to complex- valued soft SuSY-breaking parameters in the neutralino 
and chargino sector; the same in the Higgs sector, which removes the possibility of CP- 
even or CP-odd labels; the calculation of neutralino and chargino decay rates which are 
accurate for large tan/3; and matrix element weighting of particle distributions in three- 
body decays. 

2.1.5 Strong Dynamics in Electroweak Symmetry Breaking 

The simulation of strong dynamics associated with electroweak symmetry breaking in 
Pythia is based on an effective Lagrangian for the lightest resonances of a technicolor 
(TC)-like model. In TC, the breaking of a chiral symmetry in a new, strongly interacting 
gauge theory generates the Goldstone bosons necessary for electroweak symmetry break- 
ing. Bound states of technifermions provide a QCD-like spectrum of technipions (vTtc), 
technirhos (ptc), techniomegas (cutc), etc. The mass hierarchies, however, are unlike QCD 
because of the behavior of the gauge couplings in realistic models of extended TC (ETC). 
The difficulties of ETC in explaining the top quark mass while suppressing FCNC's is 
circumvented by the addition of topcolor interactions, which provide the bulk of mt. 

In ETC models, hard mass contributions to technipion masses make decays like 
Ptc TTtcTTtc kinematically inaccessible. Instead, decays like p™ — >• tt^^Wl, for exam- 
ple, dominate, where ew denotes constituent technifermions with only electroweak quan- 
tum numbers and Wl is a longitudinal W bosons. As a result, the ew technirho and 
techniomega tend to have small total widths. 

Effective couplings are derived in the valence technifermion approximation, and the 
techniparticle decays can be calculated directly [0. Technirhos and techniomegas are 
produced through kinematic mixing with gauge bosons, leading to final states containing 
Standard Model particles and/or pseudo- Goldstone bosons (technipions). 

As an additional wrinkle, SUc(3) non-singlet states are included along with the coloron 
of topcolor assisted technicolor. In this case, colored technirhos (and the coloron) can have 
substantial total widths and enhanced couplings to bottom and top quarks. 

2.2 QCD Radiation 

The matrix-element (ME) and parton-shower (PS) approaches to higher-order QCD cor- 
rections both have their advantages and disadvantages. The former offers a systematic 
expansion in orders of as, and a powerful machinery to handle multi-parton configura- 
tions on the Born level, but loop calculations are tough and lead to messy cancellations at 
small resolution scales. Resummed matrix elements may circumvent the latter problem 
for specific quantities, but then do not provide exclusive accompanying events. Parton 
showers are based on an improved leading-log (almost next-to-leading-log) approximation, 
and so cannot be accurate for well separated partons, but they offer a simple, process- 
independent machinery that gives a smooth blending of event classes (by Sudakov form 
factors) and a sensible match to hadronization. It is therefore natural to try to combine 
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these descriptions, so that ME results are recovered for widely separated partons while 
the PS sets the sub-jet structure. 

For final-state showers in Z° — > qq, where q is assumed essentially massless, such 



solutions are the standard since long |T^, e.g. by letting the shower slightly overpopulate 
the qqg phase space and then using a Monte Carlo veto technique to reduce down to the 
ME level. 

2.2.1 Initial-State Showers 

A similar technique is now available for the description of initial-state radiation in the 
production of a single colour-singlet resonance, such as ''y* /TP /W^ [|r^. The basic idea 
is to map the kinematics between the PS and ME descriptions, and to find a correction 
factor that can be applied to hard emissions in the shower so as to bring agreement with 
the matrix-element expression. The Pythia shower kinematics definitions are based on 
as the spacelike virtuality of the parton produced in a branching and z as the factor 
by which the s of the scattering subsystem is reduced by the branching. Some simple 
algebra then shows that the two qq' gW^ emission rates disagree by a factor 

D i\ - (do-/dt)ME _ + + 2m%^s 

(d(T/dt)ps + 

which is always between 1/2 and 1. The shower can therefore be improved in two ways, 
relative to the old description. Firstly, the maximum virtuality of emissions is raised from 
Qmax ~ "^w to Qmax = ^i i-^- the showcr is allowed to populate the full phase space. 
Secondly, the emission rate for the final (which normally also is the hardest) q — qg 
emission on each side is corrected by the factor R{s, t) above, so as to bring agreement 
with the matrix-element rate in the hard-emission region. In the backwards evolution 
shower algorithm ||18|, this is the first branching considered. 



The other possible 0{as) graph is qg —>■ q'W^, where the corresponding correction 
factor is 

. (d6-/dt)ME s^ + u'^ + 2mlfi 



(d(T/dt)ps (s - m^)2 + m 



w 



which lies between 1 and 3. A probable reason for the lower shower rate here is that the 
shower does not explicitly simulate the s-channel graph qg — » q* — > q'W. The g — > qq 
branching therefore has to be preweighted by a factor of 3 in the shower, but otherwise 
the method works the same as above. Obviously, the shower will mix the two alternative 
branchings, and the correction factor for a final branching is based on the current type. 

The reweighting procedure prompts some other changes in the shower. In particular, 
u < translates into a constraint on the phase space of allowed branchings, not previously 
implemented. 

Our published comparisons with data on the p^w spectrum show quite a good agree- 
ment with this improved simulation [jl^ . A worry was that an unexpectedly large primor- 



dial k±, around 4 GeV, was required to match the data in the low-p_Lw region. However, 
at that time we had not realized that the data were not fully unsmeared. The required 
primordial k± therefore drops by about a factor of two 

The method can also be used for initial-state photon emission, e.g. in the process 
e^e^ — > 7*/Z°. There the old default Q^ax = allowed no emission at large p_L, 
p± ^ niz at LEP2. This is now corrected by the increased Q^ax = ■5, and using the R of 
eq. (H) with mw — ^ mz- 
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The above method does not address the issue of next-to-leading order corrections to 
the total W cross section, which instead can be studied with more sophisticated matching 
procedures |]2D|. Also extensions to other processes can be considered in the future. 



There are also some other changes to the initial state radiation algorithm: 

• The cut on minimum gluon energy emitted in a branching is modified by an extra 
factor roughly corresponding to the I/7 factor for the boost to the hard subprocess 
frame. Earlier, when a subsystem was strongly boosted, the minimum energy re- 
quirement became quite stringent on the low-energy incoming side, and could cut 
out much radiation. 

• The angular-ordering requirement is now based on ordering p^/p rather than p^/pLi 
i.e. replacing tan6' by sin6'. Earlier the starting value (tan6')max = 10 could actually 
be violated by some bona fide emissions for strongly boosted subsystems. 

• The value of the backwards evolution of a heavy quark like c in a proton beam is 
by force kept above m^, so as to ensure that the branching g ^ cc is not 'forgotten' 
by evolving Q"^ below Qq. Thereby the possibility of having a c in the beam remnant 
proper is eliminated [^. The procedure is not forced for a photon beam, where 
charm occurs as part of the valence flavour content. 

• For incoming /x^ (or r^) beams the kinematical variables are better selected to 
represent the differences in lepton mass, and the lepton-inside-lepton parton distri- 
butions are properly defined. 

2.2.2 Final-State Showers 



The traditional final-state shower algorithm in Pythia is based on an evolution 
in = m?, i.e. potential branchings are considered in order of decreasing mass. A 
branching a — > 6c is then characterized by and z = Ef,/Ea. For the process ^* jT? — > 
qq, the first gluon emission off both q and q are corrected to the first-order matrix elements 
for 7*/Z° — s> qqg. (The and the Sudakov form factor are omitted from the comparison, 
since the shower procedure here attempts to include higher-order effects absent in the 
first-order matrix elements.) 

This matching is well-defined for massless quarks, and was originally used unchanged 
for massive ones. A first attempt to include massive matrix elements did not compensate 
for mass effects in the shower kinematics, and therefore came to exaggerate the suppression 
of radiation off heavy quarks [B^ . Now the shower has been modified to solve this issue. 



and also improved and extended to cover better a host of different reactions 



The starting point is the calculation of processes a ^ be and a — > beg, where the ratio 

Wme{xi, X2) = — ' (4) 

cr(a be) axi ax2 

gives the process-dependent differential gluon-emission rate. Here the phase space vari- 
ables are xi = 2Ei,/ma and X2 = 2Ec/ma, expressed in the rest frame of parton a. Using 
the standard model and the minimal supersymmetric extension thereof as templates, a 
wide selection of colour and spin structures have been addressed, exemplified by Z° qq, 
t bW+, H° ^ qq, t ^ bH+, Z° ^ qq, q ^ q'W+, H° ^ qq, q ^ q'H+, x ^ QQ, 
q qx, t tx, g — > qq, q — ^ qg, and t tg. The mass ratios ri = mb/nia and 
f2 = iTT'c/^a have been kept as free parameters. When allowed, processes have been cal- 
culated for an arbitrary mixture of "parities", i.e. without or with a 75 factor, like in the 
vector/axial vector structure of 7*/Z°. All the matrix elements are encoded in the new 
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function PYMAEL(NI , XI ,X2 ,R1 ,R2 , ALPHA) , where NI distinguishes the matrix elements 
and ALPHA is related to the 75 admixture. 

In order to match to the singularity structure of the massive matrix elements, the 
evolution variable is changed from rn? to rn? — rriQ^.g^eii; i-e- is the propagator 

of a massive particle. Furthermore, the z variable of a branching needs to be redefined, 
which is achieved by reducing the three-momenta of the daughters in the rest frame of 
the mother. For the shower history b ^ bg this gives a differential probability 

where the numerator 1 + of the splitting kernel for q — qg has been replaced by a 
2 in the shower algorithm. For a process with only one radiating parton in the final 
state, such as t — bW^, the ratio VFme/W^ps,i gives the acceptance probability for an 
emission in the shower. The singularity structure exactly agrees between ME and PS, 
giving a well-behaved ratio always below unity. If both b and c can radiate, there is a 
second possible shower history that has to be considered. The matrix element is here 
split in two parts, one arbitrarily associated with 6 — > 6g branchings and the other with 
c — >• eg ones. A convenient choice is PFme,i = W^me(1 + rf — r"^ — xij/x^ and Wme,2 = 
WME{^+r2 — rl—X2)/x3, which again gives matching singularity structures in WME,i/Wps,i 
and thus a well-behaved Monte Carlo procedure. 

Also subsequent emissions of gluons off the primary particles are corrected to Wme- 
To this end, a reduced-energy system is constructed, which retains the kinematics of the 
branching under consideration but omits the gluons already emitted, so that an effective 
three-body shower state can be mapped to an (xi, X2, ri, set of variables. For light 
quarks this procedure is almost equivalent with the original one of using the simple uni- 
versal splitting kernels after the first branching. For heavy quarks it offers an improved 
modelling of mass effects also in the coUinear region. 

Some further changes have been introduced, a few minor as default and some more 
significant ones as non-default options This includes the description of coherence 

effects and as arguments, in general and more specifically for secondary heavy flavour 
production by gluon splittings. 

Further issues remain to be addressed, e.g. radiation off particles with non-negligible 
width. In general, however, the new shower should allow an improved description of gluon 
radiation in many different processes. 

2.3 Beam Remnants and Multiple Interactions 
2.3.1 Beam Remnants 

In a hadron-hadron collision, the initial-state radiation algorithm reconstructs one shower 
initiator in each beam, by backwards evolution from the hard scattering. This initiator 
only takes some fraction of the total beam energy, leaving behind a beam remnant that 
takes the rest. Since the initiator is coloured, so is the remnant. It is therefore colour- 
connected to the hard interaction, and forms part of the same fragmenting system. Often 
the remnant can be complicated, e.g. a g initiator would leave behind a uud proton- 
remnant system in a colour octet state, which can conveniently be subdivided into a 
colour triplet quark and a colour antitriplet diquark, each of which are colour-connected 
to the hard interaction. The energy sharing between these two remnant objects, and their 
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relative transverse momentum, introduces additional nonperturbative degrees of freedom. 
Some of the default values have recently been updated [pT| . 

One would expect an ep event to have only one beam remnant, and an e+e" event 
none. This is not always correct, e.g. a 77 — qq interaction in an e"'"e~ event would leave 
behind the e+ and e~ as beam remnants. The photons may in their turn leave behind 
remnants. 

It is customary to assign a primordial transverse momentum to the shower initiator, 
to take into account the motion of quarks inside the original hadron, basically as required 
by the uncertainty principle. A number of the order of {k±) ^ rrip/S ~ 300 MeV could 
therefore be expected. However, in hadronic collisions much higher numbers than that 



are often required to describe data, typically of the order of 1 GeV pj, |I9[ if a Gaussian 
parameterization is used. (This number is now the default.) Thus, an interpretation as a 
purely nonperturbative motion inside a hadron is difficult to maintain. 

Instead a likely culprit is the initial-state shower algorithm. This is set up to cover 
the region of hard emissions, but may miss out on some of the softer activity, which 
inherently borders on nonperturbative physics. By default, the shower does not evolve 
down to scales below Qq = 1 GeV. Any shortfall in shower activity around or below this 
cutoff then has to be compensated by the primordial k± source, which thereby largely 
loses its original meaning. 

2.3.2 Multiple Interactions 

Multiple parton-parton interactions is the concept that, based on the composite nature 
of hadrons, several parton pairs may interact in a typical hadron-hadron collision . 



Over the years, evidence for this mechanism has accumulated, such as the recent direct 
observation by CDF |2^. The occurrences with two parton pairs at reasonably large p± 



just form the top of the iceberg, however. In the Pythia model, most interactions are at 
lower p±, where they are not visible as separate jets but only contribute to the underlying 
event structure. As such, they are at the origin of a number of key features, like the broad 
multiphcity distributions, the significant forward-backward multiplicity correlations, and 
the pedestal effect under jets. 

Since the perturbative jet cross section is divergent for p± 0, it is necessary to 
regularize it, e.g. by a cut-off at some p_Lmin scale. That such a regularization should 
occur is clear from the fact that the incoming hadrons are colour singlets — unlike the 
coloured partons assumed in the divergent perturbative calculations — and that therefore 
the colour charges should screen each other in the p± ^ limit. Also other damping 



mechanisms are possible |23]. Fits to data typically give p_Lmin ~ 2 GeV, which then 
should be interpreted as the inverse of some colour screening length in the hadron. 

One key question is the energy- dependence of p_Lmin; this may be relevant e.g. for com- 
parisons of jet rates at different Tevatron energies, and even more for any extrapolation to 
LHC energies. The problem actually is more pressing now than at the time of the original 



study |]23], since nowadays parton distributions are known to be rising more steeply at 
small X than the flat xf{x) behaviour normally assumed for small before HERA. This 
translates into a more dramatic energy dependence of the multiple-interactions rate for a 

flxed p_Lmin- 

The larger number of partons also should increase the amount of screening, however, 
as conflrmed by toy simulations |23]. As a simple flrst approximation, p±uiin is assumed 



to increase in the same way as the total cross section, i.e. with some power e ~ 0.08 [2£] 
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that, via reggeon phenomenology, should relate to the behaviour of parton distributions 
at small x and Q^. Thus the new default in PYTHIA is 



P±mUs) = (1.9 GeV) (p^j • (6) 
2.4 Fragmentation and Decays 

QCD perturbation theory, formulated in terms of quarks and gluons, is valid at short 
distances. At long distances, QCD becomes strongly interacting and perturbation theory 
breaks down. In this confinement regime, the coloured partons are transformed into 
colourless hadrons, a process called either hadronization or fragmentation. 

The fragmentation process has yet to be understood from first principles, starting 
from the QCD Lagrangian. This has left the way clear for the development of a number 
of different phenomenological models. Pythia is intimately connected with string frag- 
mentation, in the form of the time-honoured 'Lund model' [^|. This is the default for all 



applications. Improvements have been made in some areas, however. 
2.4.1 Low-Mass Strings 

A hadronic event is conventionally subdivided into sets of partons that form separate 
colour singlets. These sets are represented by strings, that e.g. stretch from a quark end 
via a number of intermediate gluons to an antiquark end. Three string mass regions may 
be distinguished for the hadronization. 

1. Normal string fragmentation. In the ideal situation, each string has a large invariant 



mass. Then the standard iterative fragmentation scheme ^ works well. In 
practice, this approach can be used for all strings above some cut-off mass of a few 
GeV. 

2. Cluster decay. If a string is produced with a small invariant mass, maybe only 
two-body final states are kinematically accessible. The traditional iterative Lund 
scheme is then not applicable. We call such a low-mass string a cluster, and consider 
it separately from above. In recent program versions, the modeling has now been 
improved to give a smooth match on to the standard string scheme in the high- 
cluster-mass limit pi|| . 

3. Cluster collapse. This is the extreme case of the above situation, where the string 
mass is so small that the cluster cannot decay into two hadrons. It is then assumed 
to collapse directly into a single hadron, which inherits the flavour content of the 
string endpoints. The original continuum of string/cluster masses is replaced by a 
discrete set of hadron masses. Energy and momentum then cannot be conserved 
inside the cluster, but must be exchanged with the rest of the event. This description 
has also been improved [pl[] . 

String systems below a threshold mass are handled by the cluster machinery. In it, an 
attempt is first made to produce two hadrons, by having the string break in the middle 
by the production of a new qq pair, with flavours and hadron spins selected according 
to the normal string rules. If the sum of the hadron masses is larger than the cluster 
mass, repeated attempts can be made to find allowed hadrons; the default is two tries. If 
an allowed set is found, the angular distribution of the decay products in the cluster rest 
framed is picked isotropically near the threshold, but then gradually more elongated along 
the string direction, to provide a smooth match to the string description at larger masses. 
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This also includes a forward-backward asymmetry, so that each hadron is preferentially 
in the same hemisphere as the respective original quark it inherits. 

If the attempts to find two hadrons fail, one single hadron is formed from the given 
flavour content. The basic strategy thereafter is to exchange some minimal amount of 
energy and momentum between the collapsing cluster and other string pieces in the neigh- 
bourhood. The momentum transfer can be in either direction, depending on whether the 
hadron is lighter or heavier than the cluster it comes from. When lighter, the excess mo- 
mentum is split off and put as an extra 'gluon' on the nearest string piece, where 'nearest' 
is defined by a space-time history-based distance measure. When the hadron is heavier, 
momentum is instead borrowed from the endpoints of the nearest string piece. 

The free parameters of the model can be tuned to data, especially to the significant 
asymmetries observed between the production of D and D mesons in vr^p collisions, with 
hadrons that share some of the tt^ flavour content very much favoured at large xp in the 



7r~ fragmentation region These spectra and asymmetries are closely related to the 
cluster collapse mechanism, and also to other effects of the colour topology of the event 
('beam drag') |2ll]. Also other parameters enter the description, however, such as the 



effective charm mass and the beam remnant structure. 
2.4.2 Baryon Production 

A new advanced scheme has been introduced for baryon production with the popcorn 
mechanism p3|, plus some minor changes to the older popcorn scheme ||3^. These new 



features currently only appear as options, with the default unchanged, and can be sepa- 
rated into three parts. 

Firstly, an improved implementation of SU(6) weights for baryon production. This 
should not be regarded as a new model, rather a more correct implementation of the old. 
However, in order to enable the user to see the effects of the SU(6) weighting separately, 
both procedures are available as different options. The main change is that, if a step 
q — i> B-|-qq' is SU(6)-rejected, the new try may now instead give a q — M-|-q' step (where 
B stands for baryon, M for meson). The old procedure leads to a slightly faster algorithm 
and a better interpretation of the input parameter for the diquark-to-quark production 
rate. However, the probability that a quark will produce a baryon and a antidiquark 
is then flavour independent, which is not in agreement with the model. Further, for 
qq — »• M -|- qq', SU(6) symmetry is included in the weights for qq', while qq is kept with 
unit probability. The procedures for qq — >^ B + q' and a final joining qq + q — >^ B are 
unchanged. 

Secondly, a suppression of diquark vertices occuring at small proper times. This is 
based on a study of the production dynamics of the three quarks that form a baryon. The 
main experimental consequence is a suppression of the baryon production rate at large 
momentum fraction. This in particular implies a smaller rate of first-rank light baryon 
production, while charm and bottom baryons are less affected (since the production proper 
time is larger for a heavy hadron than a light one of the same momentum). It thereby 
substitutes and explains the older brute-force possibility to suppress the production of 
first-rank baryons. 

Thirdly, a completely new fiavour algorithm for baryons and popcorn mesons, also 
using the small-proper-time suppression above. While the old popcorn alternative allowed 
at most one meson to be produced in between the baryon and the antibaryon, the new 
model allows an arbitrary number. The new fiavour model makes explicit use of the 
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popcorn suppression factor exp(— 2m_|_M_|_/K), where m_L is the transverse mass of the 
quark creating the colour fluctuation, M_|_ is the total invariant transverse mass of the 
popcorn meson system, and k is the string tension constant. Thus two parameters, 
representing the mean 2m^/ n for light quarks and s-quarks, respectively, govern both 
diquark and popcorn meson production. A corresponding parameter is introduced for the 
fragmentation of strings that contain diquarks already from the beginning, i.e. baryon 
remnants. The new procedure therefore requires far fewer parameters than the old one, 
and still provides a comparable quality in the description of the various baryon production 
rates. This was investigated in detail in (The concluding worry of an "improper 



treatment" was caused by an unfortunate misunderstanding and can be disregarded.) 
Other features, such as baryon correlations, are also modified. 

Several new routines have been added, and the diquark code has been extended with 
information about the curtain quark flavour, i.e. the qq pair that is shared between the 
baryon and antibaryon, but this is not visible externally. Some parameters are no longer 
used, while others have to be given modified values, as described in the long writeup. 

2.4.3 Interconnection Effects 

The widths of the W, Z and t are all of the order of 2 GeV. A standard model Higgs with 
a mass above 200 GeV, as well as many supersymmetric and other beyond the standard 
model particles would also have widths in the multi-GeV range. Not far from threshold, 
the typical decay times r = l/F ^ 0.1 fm <^ Thad ~ 1 fm. Thus hadronic decay systems 
overlap, between a resonance and the underlying event, or between pairs of resonances, 
so that the final state may not contain independent resonance decays. 

So far, studies have mainly been performed in the context of W pair production at 
LEP2. Pragmatically, one may here distinguish three main eras for such interconnection: 

1. Perturbative: this is suppressed for gluon energies uo > T hy propagator/timescale 
effects; thus only soft gluons may contribute appreciably. 

2. Nonperturbative in the hadroformation process: normally modeled by a colour re- 
arrangement between the partons produced in the two resonance decays and in the 
subsequent parton showers. 

3. Nonperturbative in the purely hadronic phase: best exemplified by Bose-Einstein 
effects. 

The above topics are deeply related to the unsolved problems of strong interactions: 
confinement dynamics, I/Nq effects, quantum mechanical interferences, etc. Thus they 
offer an opportunity to study the dynamics of unstable particles, and new ways to probe 
confinement dynamics in space and time p6| , |37| , hut they also risk to limit or even spoil 
precision measurements. 

The reconnection scenarios outlined in are now available, plus also an option 



along the lines suggested in |3g]. Currently they can only be invoked in process 25, 
e"'"e~ W"'"W~ qiq2l3Q45 which is the most interesting one for the foreseeable future. 
(Process 22, e+e~ — >• 7*/Z° 7*/Z° (li^i'^i^i can also be used, but the travel distance 
is calculated based only on the TP propagator part.) If normally the event is considered 
as consisting of two separate colour singlets, qiq2 from the W"^ and q3q4 from the W~, a 
colour rearrangement can give two new colour singlets qiq4 and q3q2. It therefore leads 
to a different hadronic final state, although differences usually turn out to be subtle and 
difficult to isolate [^. When also gluon emission is considered, the number of potential 
reconnection topologies increases. Apart from the overall rate of reconnection, the sce- 
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narios in Pythia differ in the relative probability assigned to each of these topologies, 
based on their properties in momentum space and/or space-time. For instance, scenario 
I is based on an analogy with type I superconductors, with the colour field represented 
by extended flux tubes. By contrast, scenario II assumes that narrow vortex lines carry 
all the topological information, like in type II superconductors, even if the full energy is 
stored over a wider region. 

Bose-Einstein effects are simulated in a simplified manner, by introducing small mo- 
mentum shifts in identical final-state mesons (primarily and 7r°) so as to bring them 
closer to each other The shifts can be chosen to reproduce a desired BE enhancement 



shape for small relative momentum Q = ymfj — 4mf between identical bosons i and j. 
Typically the shape is chosen as a Gaussian, f2{Q) = 1 -l- A exp(— Q^i?^), with A and R 
two free parameters. The input is only exactly reproduced in the limit of an isotropic and 
low-density initial particle distribution; since these conditions are not completely fulfilled 
in reality, there are distortions for better or worse. (The nontrivial three-particle 



correlations in data are described qualitatively, although not quantitatively.) 

A major shortcoming of the algorithm is that energy is not automatically conserved, 
even though three-momentum is. In the original algorithm, this was solved by a uniform 
rescaling of all three-momenta, with undesirable side effects e.g. when studying BE effects 
in W"'"W~ hadronic final states. In the current version, several new options have been 
added that, based on different principles, instead shifts pairs apart. The default one, 
BE32, operates on identical particles, introducing an extra factor 

1 + aXexpi-Q^Ryg) {l - exp{-Q'^RyA)} (7) 

to /2(<5)- Here a is a negative number adjusted event by event for overall energy con- 
servation, with (a) ~ —0.25. This scenario can be viewed as a simplified version of a 
dampened oscillating correlation function, where only the first peak and dip has been re- 
tained. Further new options have also been introduced specifically geared towards studies 
of W"'"W~ hadronic events, e.g. to include the effects of the separated W"*" and W~ decay 
vertices. 



2.4.4 Decays 

Two separate decay treatments exist in Pythia. One is making use of a set of tables 
where branching ratios and decay modes are stored, and is used e.g. for hadronic decays, 
where branching ratios normally cannot be calculated from first principles. 

The other treatment is used for a set of fundamental resonances in or beyond the 
standard model, such as t, Z°, W^, h*^, super symmetric particles, and many more. Char- 
acteristic here is that these resonances have perturbatively calculable widths to each of 
their decay channels. The decay products are typically quarks, leptons, or other reso- 
nances. In decays to quarks, parton showers are automatically added to give a more 
realistic multijet structure, and one may also allow photon emission off leptons. If the 
decay products in turn are resonances, further decays are necessary. Often spin informa- 
tion is available in resonance decay matrix elements, leading to nonisotropic decays. This 
part has been improved in several processes, but is still missing in many others. 

The routine used to calculate the partial and total width of resonances (now expressed 
in GeV throughout), has been expanded for all the new particles and decay modes in- 
troduced. Some alternative calculation schemes have also been adopted, e.g. based on a 



18 



simple rescaling of the on-shell widths rather than a complete recalculation (which may 
at times not be feasible) based on the current mass. 

The width to be used in the denominator of a resonance propagator is only well- 
defined near the peak. Well away from the peak, an unfortunate choice may lead to a loss 
of cancellation between resonant and nonresonant diagrams. A special problem exists for 
a massive standard model Higgs, where the width Fh oc is so large that the choice of 
s dependence of the width significantly influences the resonance peak shape. Following 
H^, the default now is Fh oc m^-\/I. 



3 Program News 

Essentially all of the basic philosophy and framework remain from the previous Pythia 
and Jetset versions, so no user familiar with these should feel at loss with Pythia 6.1. 
Most of the changes and additions instead are under the surface, and are only visible as 
new options added to the existing repertoire. However, some changes are fairly obvious, 
and other less obvious ones still of general interest. These will be covered in this section, in 
fairly general terms. Again we refer to the Pythia web page for a detailed documentation. 



3.1 Coding conventions 

As before, the Fortran 77 standard is adhered to. A very few minor extensions may be 



used in isolated places, like the 7-character names of the Pdflib routines W^, but are 
not known to cause problems on any compiler in use. 

An obvious consequence of the Pythia/Jetset code merging is that the old Jetset 
routines and commonblocks have been renamed to begin with PY (instead of LU or UL), 
just like the Pythia ones. In most cases, the rest of the name is unchanged, but there 
are a few exceptions, mainly RLU^PYR, KLU^PYK, PLU^PYP and LUXTOT^PYXTEE. Three 
integer functions now begin with PY, namely PYK, PYCHGE and PYCOMP, and therefore have 
to be declared extra. The LUDATA block data has been merged into PYDATA, and the test 
routine LUTEST into PYTEST. For rotations and boosts, the PYROBO routine now requires 
the range of affected entries to be given, like the old LUDBRB but unlike LUROBO (but 0,0 
as range arguments gives back the old LUROBO behaviour). 

All real variables are now in DOUBLE PRECISION, which is assumed to mean 64 bits, 
and also real constants have been promoted to the higher precision. This is required to 
ensure proper functioning at currently studied energies, such as the LHC and beyond. To 
take into account this, all routines begin with the declarations 

C... Double precision and integer declarations. 
IMPLICIT DOUBLE PRECISION (A-H, 0-Z) 
IMPLICIT INTEGER(I-N) 
INTEGER PYK, PYCHGE, PYCOMP 

and users should do the same in their main programs. 

On a machine where DOUBLE PRECISION would give 128 bits, it may make sense to 
use compiler options to revert to 64 bits, since the program is anyway not constructed to 
make use of 128 bit precision. 

The random number generator is the same as in previous versions |^3|, but has now 
been expanded to operate with a 48 bit mantissa for the real numbers. 
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Table 4: The current form of the main commonblock declarations. 

COMMON/PYJETS/N , NPAD , K (4000 , 5) , P (4000 , 5) , V (4000 , 5) 
C0MM0N/PYDAT1/MSTU(200) ,PARU(200) ,MSTJ(200) ,PARJ(200) 
C0MM0N/PYDAT2/KCHG(500,4) ,PMAS(500,4) ,PARF(2000) ,VCKM(4,4) 
C0MM0N/PYDAT3/MDCY(500,3) ,MDME(4000 , 2) , BRAT (4000) ,KFDP(4000 , 5) 
C0MM0N/PYDAT4/CHAF (500 , 2) 
CHARACTER CHAF*16 
C0MM0N/PYDATR/MRPY(6) ,RRPY(100) 

COMMON/PYSUBS/MSEL , MSELPD , MSUB (500) , KFIN (2 , -40 : 40) , CKIN (200) 
C0MM0N/PYPARS/MSTP(200) ,PARP(200) ,MSTI(200) , PARI (200) 
C0MM0N/PYINT1/MINT(400) ,VINT(400) 

C0MM0N/PYINT2/ISET(500) , KFPR(500 , 2) , COEF (500 , 20) , IC0L(40 , 4, 2) 
C0MM0N/PYINT3/XSFX(2,-40:40) ,ISIG(1000,3) ,SIGH(1000) 
C0MM0N/PYINT4/MWID (500) , WIDS (500 , 5) 
C0MM0N/PYINT5/NGENPD , NGEN (0 : 500 , 3) , XSEC (0 : 500 , 3) 
C0MM0N/PYINT6/PR0C (0 : 500) 
CHARACTER PR0C*28 

COMMON/PYMSSM/IMSS (0 : 99) , RMSS (0 : 99) 

COMMON/PYUPPR/NUP , KUP (20 , 7) , NFUP , IFUP ( 10 , 2) , PUP (20 , 5) , Q2UP (0:10) 
C0MM0N/PYBINS/IHIST(4) ,INDX(1000) ,BIN(20000) 



Fortran 77 makes no provision for double-precision complex numbers, but since 
COMPLEX is used only sparingly, no problems should be expected from this omission. For 
the technicolor processes, some variables are declared COMPLEX* 16 in the PYSIGH routine. 
Should the compiler not accept this, that one declaration can be changed to COMPLEX with 
some drop in precision for the affected processes. 

Several compilers report problems when an odd number of integers precede a double- 
precision variable in a commonblock. Therefore an extra integer has been introduced as 
padding in a few instances (NPAD, MSELPD and NGENPD in Table |). 

In order to cater for the increased offering of subprocesses, some arrays in common- 
blocks have been expanded. A few, such as PYINT4, have also been reorganized to represent 
improvements in the physics modeling. Most commonblocks and commonblock variables 
are easily recognizable from previous program versions, however. The current complement 
is given in Table ^, omitting some of the less interesting ones. 

Since Fortran 77 provides no date-and-time routine, PYTIME allows a system-specific 
routine to be interfaced, with some commented-out examples given in the code. This 
routine is only used for cosmetic improvements of the output, however, so can be left at 
the default with time given. 

For a program written to run Pythia 5 and Jetset 7, most of the conversion required 
for Pythia 6 is fairly straightforward, and can be automatized. Both a simple Fortran 
routine and a more sophisticated Perl |44] script exists to this end. Some manual checks 
and interventions may still be required. 
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Table 5: New or modified particle codes or names. 



Renamed: 




7 


b' 


8 


t' 


17 


t' 


18 


K 


25 




35 


HO 


Moved: 




100443 


^' 


100553 




4000001 


d* 


4000002 


u* 


4000011 


e* 


4000012 


< 




SuSY: 




1000001 


di 


1000002 




1000003 


sl 


1000004 


Cl 


1000005 


bi 


1000006 


ti 


1000011 


Gl 


1000012 




1000013 


Ur, 


1000014 




1000015 


fi 


1000016 




1000021 


g 


1000022 




1000023 




1000024 


xt 



Susy: 




2000001 


di? 


2000002 




2000003 


Si? 


2000004 




2000005 


b2 


2000006 


t2 


2000011 


eji 


2000012 




2000013 


Ur 


2000014 




2000015 


T2 


2000016 


VrR 


1000025 


xl 


1000035 


xl 


1000037 


xt 


1000039 


G 



3.2 Particle codes and data 

A number of new particle codes KF have been introduced, or modified, see Table |^. Mostly 



this is based on the PDG-agreed conventions [45, 46 , but some not yet standardized codes 



appear in the 'empty' range 41-80. Furthermore, the fourth generation fermions and 
neutral scalar Higgs states have been renamed. The two fermion spartners are labelled 
left and right, except in the third generation, where an expected larger mixing makes the 
two mass eigenstates a better choice of classification. 

The top hadrons are gone. It is now known that top is too short-lived to form hadronic 
bound states, so a reasonable description is instead to have the top quarks decay before 
hadronization is considered. The same is now assumed about a hypothetical fourth gen- 
eration. Should the need ever arise in the future to consider a new long-lived coloured 
object, an effective description of a hadron as a small string with an ordinary colour- 
matching flavour at the other end should be sufficient. One such example would be 



leptoquark-hadrons ^ 



Bottom hadrons are now defined individually, e.g. the previous common decay scheme 
is gone in favour of individual branching ratios for each hadron. On the other hand, given 
the sketchy knowledge of many branching ratios, the default description is still fairly 
standardized. 

Decay data is mainly based on the 1996 PDG edition but with many 'educated 
guesses' to fill in missing information. 

Since running fermion masses are used in an increasing number of processes, e.g. for 
Higgs couplings, a function PYMRUN(KF,Q2) has been introduced to give the mass as a 
function of scale. 
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The compressed codes, KC = PYCOMP(KF), are completely changed. We remind that 
KF can range up to seven-digit codes, plus a sign. They therefore cannot be used to 
directly access information in particle data tables. The KC codes range between 1 and 
500, and give the index to the particle data arrays. Each KF code is now one-to-one 
associated with a KC code; the only ambiguity is that KC does not distinguish antiparticles 
from particles. Whereas KF codes below 100 still obey KF = KC, the mapping of codes 
above 100 is completely changed. It is no longer hard-coded in PYCOMP, but defined by 
the fourth component of the KCHG array. Therefore it can be changed or expanded during 
the course of a run, either by PYUPDA calls or by direct user intervention. 

3.3 New Options 

A large amount of new options have been added, related to almost all the physics changes 
above and more, and we here only mention some of the more significant ones. 

The inclusion of SuSY processes means that all the SPythia PYMSSM commonblock 
switches and parameters are inherited. New parameters are added also for other new 
physics scenarios, such as technicolor and doubly-charged Higgses. 

The extensions to the physics of virtual photons, outlined in section |2.1.3|, has resulted 
in two sets of new possibilities. One is in the description of the virtual-photon flux, where 
new CKIN switches has been introduced, e.g. to set the range of photon x and Q^. This 
is available when PYINIT is called with ' gamma/ lept on' as beam or target, to denote 
that the photon flux inside the lepton has to be considered as a new administrative layer, 
also documented in the event record. The other is the new physics machinery. Here the 
main switch is MSTP(14) that sets the assumed nature of the photon or photons, e.g. 
'a direct photon from the left collides with a VMD one from the right'. The default is 
the most general mixture, meaning 4 components for 7*p and 13 for 7*7*. This is the 
relevant approach for studies of QCD processes. There is no corresponding automatic 
mixing machinery for other processes, so then the relevant contributing components have 
to be handled separately and added afterwards. Further options are available for several 
of the components, e.g. the DIS process dampening in the — > limit, the relative 
normalization of the GVMD spectrum, the scale choice for parton distributions, and the 
possibility to add the effects of a longitudinal resolved contribution. 

The matrix-element options for e+e^ 7*/Z° ^ 2, 3 or 4 partons have previously 
only been available via the LUEEVT/PYEEVT routine, that suffers from problems of its own 
in having a rather old-fashioned machinery for QED initial-state radiation and electroweak 
parameters. Now the QCD matrix-element description is accessible as an option to the 
shower default for e"'"e~ events generated with subprocess 1 of the standard Pythia 
machinery. 

3.4 Interfaces 

While Pythia contains an extensive library of subprocesses, it is far from up to all the 
requirements of the experimental community. Both further processes and a more detailed 
treatment of the existing ones is required at times. In particular, it is not uncommon 
with a generator dedicated to one specific process, where also higher-order electroweak 
corrections, absent in Pythia, have been included in the cross section. None of these 
programs are geared to handle the QCD aspects of parton showers and hadronization, 
however, so it makes sense to combine the individual strengths. 
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A generic facility to include external processes exists since long in Pythia. Here one 
can feed in partonic configurations from an external generator, together with some basic 
information on colour flow and which partons are allowed to radiate, and let Pythia 
construct a complete event based on this information. For the simple configurations 
encountered in e^e~ annihilation events, this would often be overkill, since neither the 
initial-state QCD radiation nor beam-remnant treatment of the generic (hadronic) colli- 
sion is present. 

Based on the concepts presented in the LEP2 workshop |^ , a few simpler alternatives 
are therefore now provided for this kind of tasks: 

• CALL PY2FRM(IRAD,ITAU,IC0M) allows a parton shower to develop and partons to 
hadronize from a given two-fermion starting point. IRAD sets whether quarks are 
allowed also to radiate photons or not, ITAU whether r leptons should be decayed or 
not, and ICOM whether the input and output event record is HEPEVT or PYJETS. An 
arbitrary number of photons (e.g. from initial-state radiation) may also be stored 
with the input. 

• CALL PY4FRM(AT0TSQ,A1SQ,A2SQ,ISTRAT, IRAD, ITAU, ICOM) allows parton show- 
ers to develop and partons to hadronize from a given four-fermion starting point. 
The extra parameters can be used to select between the two colour pairings allowed 
for a qiq2q3q4 state, according to some different strategies when interference terms 
do not allow unique probabilities to be found. 

• CALL PY6FRM (P 12, P 13, P21,P23,P31,P32,PT0P, IRAD, ITAU, ICOM) allows parton 
showers to develop and partons to hadronize from a given six-fermion starting point. 
The Pi j parameters give the relative probabilities for the six colour pairings allowed 
for a six-quark state, and PTOP the probability that the event originates from a tt 
pair (in which case the shower handling has to be different than e.g. in a Z°W"'"W~ 
event). 

The above routines are not set up to handle QCD four-jet events, i.e. events of 
the types qqgg and qqq'q', with q'q' coming from a gluon branching. Such events are 
generated in normal parton showers, but not necessarily at the right rate (a problem 
that may be especially interesting for massive quarks like b). Therefore one would like 
to start a QCD parton shower from a given four-parton configuration. Some time ago, a 



machinery was developed to handle this kind of occurences [49 . This approach has now 



been adapted to the current Pythia version, in a somewhat modified form. In it, an 
imagined shower history of two branchings is (re) constructed from the four-parton state, 
according to relative probabilities derived in the shower language. Thereafter a normal 
shower is allowed to develop, with branchings chosen at random except for these two 
predetermined ones. The routine CALL PY4 JET (PMAX, IRAD, ICOM) takes an original four- 
parton configuration stored in HEPEVT or PYJETS and lets a shower develop as described 
above. PMAX can be used to set the maximum virtuality of those parts of the shower not 
given from the parton configuration itself, either to a fixed value or to the lowest virtuality 
of the reconstructed shower. 



3.5 Utilities 

The clustering algorithm PYCLUS has been extended also to accept the Durham distance 
30[] as an alternative. This is p^-based, like the original LUCLUS distance measure. 



measure 



but differs in the details. 

The GBOOK histogramming package was written in 1979 as a lightweight substitute 
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for HBOOK [^] before that program was available in Fortran 77. The one-dimensional 
histogram part now appears in the standard distribution, in order to make the sample 
runs offered on the web a bit more realistic. The main routines are: 

• CALL PYBOOK ( ID , TITLE , NX , XL , XU) to book a one-dimensional histogram with inte- 
ger identifier ID (in the range 1 — 1000), character title TITLE and NX bins stretching 
from XL to XU. 

• CALL PYFILL(ID,X,W) to fill histogram ID at position X with weight W. 

• CALL PYFACT(ID,F) to rescale the contents of histogram ID by a factor F. 

• CALL PY0PER(ID1,0PER,ID2,ID3,F1,F2) to perform operations on several his- 
tograms, such as adding or dividing them by each other. 

• CALL PYDUMP(MDUMP,LFN,NHI,IHI) to dump histogram contents to a file from 
which they could be read in for plotting in another program. 

• CALL PYHIST to print all histograms in a simple line-printer mode, and thereafter 
reset histogram contents. 

A commonblock of dimension 20000 is used to store the histograms; this size may need 
to be expanded if many histograms are to be booked. 

The PYUPDA routine has been expanded with a new option that allows a set of particle 
data to be read in, in tabular form as before, as an addition to or partial replacement of 
the existing particle data. 



4 Summary and Outlook 

We have here given a very brief survey of news in the Pythia 6.1 program. A more 
detailed description of physics and programs is available separately . Any serious user 
should turn to this publication, and to the original physics papers, for further information. 

The treasure trove for information is the Pythia webpage, 

http : //www. thep . lu. se/~torbjorn/Pythia.html , 
where one may find the current and previous subversions, with documentation, sample 
main programs, links to related programs, etc. 

The Pythia program is continuously being developed. We are aware of many physics 
shortcomings, which hopefully will be addressed in the future. It is in the nature of a 
program of this kind never to be finished, at least as long as it is of importance for the 
high-energy physics experimental community. 

The main visible change in the future is the transition to C++ as the programming 
language for Pythia 7. Even if much of the physics will be carried over unchanged, 
none of the existing code will survive. The structure of the event record and the whole 
administrative apparatus is completely different from the current one, in order to allow 
a much more general and flexible formulation of the event generation process. Following 
the formulation of a strategy document ||5^, a first proof-of-concept version was released 



recently |£4|. So far it only contains one reasonably complete physics module, however, 
namely that of string fragmentation. More realistic versions should follow, but it will 
take a long time to convert all important physics components from Pythia 6. The two 
versions therefore will coexist for several years, with the Fortran one used for physics 
'production' and the C++ one for exploration of the object-oriented approach that will 
be standard at the LHC. 
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